
#ifndef _BILINEAR_H
#define _BILINEAR_H_

#include "Tuna.hpp"

class Bilinear 
{
///
///  Traits to define the type of arrays to be used.
///
///  floatTinyArray_t : Array of floats mostly used for deltas of meshes
///    
    typedef TunaArray<double, 2>::tiny floatTinyArray_t;
    typedef TunaArray<double, 2>::huge ScalarField2D;
    
 public:
    Bilinear() { };
    ~Bilinear() { };

/*!
 * This function locate the position of the particle on the 2D mesh. 
 * The cell where the particle is located is identified with the indexes 
 * (I, J). The local position inside the cell (I,J) is defined by 
 * (alpha, beta).
 * \param pos Position of the particle.
 * \param deltas Deltas on each axis of the mesh.
 */
    inline void locate(floatTinyArray_t &pos, const floatTinyArray_t &deltas)
	{
	    double a = pos(0) / deltas (0);
	    double b = pos(1) / deltas (1);
	    
	    I = static_cast<int>(a);
	    J = static_cast<int>(b);
	    
	    alpha = a - I;
	    beta  = b - J;
	}

/*!
 * This is a bilinear interpolation of the velocity to the actual position
 * of the particle. The result is stored in the array UVW.
 * \param pos Position of the particle.
 * \param deltas Deltas on each axis of the mesh.
 * \param u X-component of the velocity.
 * \param v Y-component of the velocity.
 */   
    inline floatTinyArray_t interpolate(floatTinyArray_t &pos,
					const floatTinyArray_t &deltas,
					ScalarField2D &u,
					ScalarField2D &v)
	{
	    locate(pos, deltas);
	    
	    UVW(0) = ( u(I  ,J  ) * (1.0 - alpha) +
		       u(I+1,J  ) * alpha ) * (1.0 - beta) +
		     ( u(I  ,J+1) * (1.0 - alpha) +
		       u(I+1,J+1) * alpha ) * beta;
	    
	    UVW(1) = ( v(I  ,J  ) * (1.0 - alpha) +
		       v(I+1,J  ) * alpha ) * (1.0 - beta) +
		     ( v(I  ,J+1) * (1.0 - alpha) +
		       v(I+1,J+1) * alpha ) * beta;
	    return UVW;
	}

 private:
    int I, J;
    double alpha, beta;
    floatTinyArray_t UVW;
};

#endif // _BILINEAR_H_
